function [ys,rw,res]=marray_lowess2(x,y,f,nsteps,delta)
%[ys,rw,res]=lowess(x,y,f,nsteps,delta)
% x: must be ordered from smallest to largest. 
% y: cordinates of the points on the scatterplot.
% f: is the fraction of points used to compute each fitted value
% as f increases the smoothed values become smoother, choosing f from 0.2 to
% 0.5 should be solve most of the problems
% nsteps: the number of iterations in the robust fit. nsteps=0 nonrobust fit is 
% returned. nsteps=2 should server most purpose.
% delta: nonnegative parameter which may be used to save computations; if n is less than
% 100, set delta to 0, if n greater than 100 you should try by yourself, for
% example, delta=length(x)/k, k=50.
%
n=length(y);
if n<2
   ys(1)=y(1);
end
ns=max(min(fix(f*n),n),2);
for iter=1:nsteps+1
   nleft=1;
   nright=ns;
   last=0;
   i=1;
   while last<n
   	while nright <n
      		d1=x(i)-x(nleft);
      		d2=x(nright+1)-x(i);
      		if d1>d2
        		nleft=nleft+1;
        		nright=nright+1; 
      		else
        		break;
      		end  
   	end
   
   	if iter>1
    		%iter, rw(1:5)
    		[ys(i),res,ok]=marray_lowest(x,y,n,x(i),nleft,nright,iter>1,rw);
   	else
     		[ys(i),res,ok]=marray_lowest(x,y,n,x(i),nleft,nright,iter>1);
   	end
  	if ~ok
      		ys(i)=y(i);
   	end
   	if last<i-1
      		denom=x(i)-x(last);
      		for j=last+1:i-1
         		alpha=(x(j)-x(last))/denom;
         		ys(j)=alpha*ys(i)+(1-alpha)*ys(last);
      		end
   	end
   	last=i;
   	cut=x(last)+delta;
   	for i=last+1:n
      		if x(i)>cut 
         		break;
      		end
      		if x(i)==x(last)
         		ys(i)=ys(last);
         		last=i;
      		end
   	end
   	i=max(last+1,i-1);
   end
	
   for i=1:n
   	res(i)=y(i)-ys(i);
   end

   if iter>nsteps
   	break;
   end

   for i=1:n
   	rw(i)=abs(res(i));
   end
   rw=sort(rw); %Check

   m1=n/2+1;
   m2=n-m1+1;
   m1=ceil(m1);
   m2=ceil(m2);
   cmad=3*(rw(m1)+rw(m2));
   c9=0.999*cmad;
   c1=0.001*cmad;
   for i=1:n
   	r=abs(res(i));
   	if r<=c1
      		rw(i)=1;
   	elseif r>c9
      		rw(i)=0;
   	else
      		rw(i)=(1-(r/cmad)^2)^2;
   	end
    end
end

      

   

   
         
   


